Identification of Interpretable Clusters and Associated Signatures in Breast Cancer Single-Cell Data: A Topic Modeling Approach.

Topic modeling is a popular technique in machine learning and natural language processing, where a corpus of text documents is classified into themes or topics using word frequency analysis. This approach has proven successful in various biological data analysis applications, such as predicting cancer subtypes with high accuracy and identifying genes, enhancers, and stable cell types simultaneously from sparse single-cell epigenomics data. The advantage of using a topic model is that it not only serves as a clustering algorithm, but it can also explain clustering results by providing word probability distributions over topics. Our study proposes a novel topic modeling approach for clustering single cells and detecting topics (gene signatures) in single-cell datasets that measure multiple omics simultaneously. We applied this approach to examine the transcriptional heterogeneity of luminal and triple-negative breast cancer cells using patient-derived xenograft models with acquired resistance to chemotherapy and targeted therapy. Through this approach, we identified protein-coding genes and long non-coding RNAs (lncRNAs) that group thousands of cells into biologically similar clusters, accurately distinguishing drug-sensitive and -resistant breast cancer types. In comparison to standard state-of-the-art clustering analyses, our approach offers an optimal partitioning of genes into topics and cells into clusters simultaneously, producing easily interpretable clustering outcomes. Additionally, we demonstrate that an integrative clustering approach, which combines the information from mRNAs and lncRNAs treated as disjoint omics layers, enhances the accuracy of cell classification.

PhysiBoSS 2.0: a sustainable integration of stochastic Boolean and agent-based modelling frameworks.

In systems biology, mathematical models and simulations play a crucial role in understanding complex biological systems. Different modelling frameworks are employed depending on the nature and scales of the system under study. For instance, signalling and regulatory networks can be simulated using Boolean modelling, whereas multicellular systems can be studied using agent-based modelling. Herein, we present PhysiBoSS 2.0, a hybrid agent-based modelling framework that allows simulating signalling and regulatory networks within individual cell agents. PhysiBoSS 2.0 is a redesign and reimplementation of PhysiBoSS 1.0 and was conceived as an add-on that expands the PhysiCell functionalities by enabling the simulation of intracellular cell signalling using MaBoSS while keeping a decoupled, maintainable and model-agnostic design. PhysiBoSS 2.0 also expands the set of functionalities offered to the users, including custom models and cell specifications, mechanistic submodels of substrate internalisation and detailed control over simulation parameters. Together with PhysiBoSS 2.0, we introduce PCTK, a Python package developed for handling and processing simulation outputs, and generating summary plots and 3D renders. PhysiBoSS 2.0 allows studying the interplay between the microenvironment, the signalling pathways that control cellular processes and population dynamics, suitable for modelling cancer. We show different approaches for integrating Boolean networks into multi-scale simulations using strategies to study the drug effects and synergies in models of cancer cell lines and validate them using experimental data. PhysiBoSS 2.0 is open-source and publicly available on GitHub with several repositories of accompanying interoperable tools.

Corrigendum to "Modeling signaling pathways in biology with MaBoSS: From one single cell to a dynamic population of heterogeneous interacting cells" [Comput. Struct. Biotechnol. 20 (2022) 5661-5671].

[This corrects the article DOI: 10.1016/j.csbj.2022.10.003.].

Multiscale model of the different modes of cancer cell invasion.

MOTIVATION: Mathematical models of biological processes altered in cancer are built using the knowledge of complex networks of signaling pathways, detailing the molecular regulations inside different cell types, such as tumor cells, immune and other stromal cells. If these models mainly focus on intracellular information, they often omit a description of the spatial organization among cells and their interactions, and with the tumoral microenvironment. RESULTS: We present here a model of tumor cell invasion simulated with PhysiBoSS, a multiscale framework, which combines agent-based modeling and continuous time Markov processes applied on Boolean network models. With this model, we aim to study the different modes of cell migration and to predict means to block it by considering not only spatial information obtained from the agent-based simulation but also intracellular regulation obtained from the Boolean model.Our multiscale model integrates the impact of gene mutations with the perturbation of the environmental conditions and allows the visualization of the results with 2D and 3D representations. The model successfully reproduces single and collective migration processes and is validated on published experiments on cell invasion. In silico experiments are suggested to search for possible targets that can block the more invasive tumoral phenotypes. AVAILABILITY AND IMPLEMENTATION: https://github.com/sysbio-curie/Invasion_model_PhysiBoSS.

Modeling signaling pathways in biology with MaBoSS: From one single cell to a dynamic population of heterogeneous interacting cells.

As a result of the development of experimental technologies and the accumulation of data, biological and molecular processes can be described as complex networks of signaling pathways. These networks are often directed and signed, where nodes represent entities (genes/proteins) and arrows interactions. They are translated into mathematical models by adding a dynamic layer onto them. Such mathematical models help to understand and interpret non-intuitive experimental observations and to anticipate the response to external interventions such as drug effects on phenotypes. Several frameworks for modeling signaling pathways exist. The choice of the appropriate framework is often driven by the experimental context. In this review, we present MaBoSS, a tool based on Boolean modeling using a continuous time approach, which predicts time-dependent probabilities of entities in different biological contexts. MaBoSS was initially built to model the intracellular signaling in non-interacting homogeneous cell populations. MaBoSS was then adapted to model heterogeneous cell populations (EnsembleMaBoSS) by considering families of models rather than a unique model. To account for more complex questions, MaBoSS was extended to simulate dynamical interacting populations (UPMaBoSS), with a precise spatial distribution (PhysiBoSS). To illustrate all these levels of description, we show how each of these tools can be used with a running example of a simple model of cell fate decisions. Finally, we present practical applications to cancer biology and studies of the immune response.

BIODICA: a computational environment for Independent Component Analysis of omics data.

SUMMARY: We developed BIODICA, an integrated computational environment for application of independent component analysis (ICA) to bulk and single-cell molecular profiles, interpretation of the results in terms of biological functions and correlation with metadata. The computational core is the novel Python package stabilized-ica which provides interface to several ICA algorithms, a stabilization procedure, meta-analysis and component interpretation tools. BIODICA is equipped with a user-friendly graphical user interface, allowing non-experienced users to perform the ICA-based omics data analysis. The results are provided in interactive ways, thus facilitating communication with biology experts. AVAILABILITY AND IMPLEMENTATION: BIODICA is implemented in Java, Python and JavaScript. The source code is freely available on GitHub under the MIT and the GNU LGPL licenses. BIODICA is supported on all major operating systems. URL: https://sysbio-curie.github.io/biodica-environment/.

UPMaBoSS: A Novel Framework for Dynamic Cell Population Modeling.

Mathematical modeling aims at understanding the effects of biological perturbations, suggesting ways to intervene and to reestablish proper cell functioning in diseases such as cancer or in autoimmune disorders. This is a difficult task for obvious reasons: the level of details needed to describe the intra-cellular processes involved, the numerous interactions between cells and cell types, and the complex dynamical properties of such populations where cells die, divide and interact constantly, to cite a few. Another important difficulty comes from the spatial distribution of these cells, their diffusion and motility. All of these aspects cannot be easily resolved in a unique mathematical model or with a unique formalism. To cope with some of these issues, we introduce here a novel framework, UPMaBoSS (for Update Population MaBoSS), dedicated to modeling dynamic populations of interacting cells. We rely on the preexisting tool MaBoSS, which enables probabilistic simulations of cellular networks. A novel software layer is added to account for cell interactions and population dynamics, but without considering the spatial dimension. This modeling approach can be seen as an intermediate step towards more complex spatial descriptions. We illustrate our methodology by means of a case study dealing with TNF-induced cell death. Interestingly, the simulation of cell population dynamics with UPMaBoSS reveals a mechanism of resistance triggered by TNF treatment. Relatively easy to encode, UPMaBoSS simulations require only moderate computational power and execution time. To ease the reproduction of simulations, we provide several Jupyter notebooks that can be accessed within the CoLoMoTo Docker image, which contains all software and models used for this study.

Patient-specific Boolean models of signalling networks guide personalised treatments.

Prostate cancer is the second most occurring cancer in men worldwide. To better understand the mechanisms of tumorigenesis and possible treatment responses, we developed a mathematical model of prostate cancer which considers the major signalling pathways known to be deregulated. We personalised this Boolean model to molecular data to reflect the heterogeneity and specific response to perturbations of cancer patients. A total of 488 prostate samples were used to build patient-specific models and compared to available clinical data. Additionally, eight prostate cell line-specific models were built to validate our approach with dose-response data of several drugs. The effects of single and combined drugs were tested in these models under different growth conditions. We identified 15 actionable points of interventions in one cell line-specific model whose inactivation hinders tumorigenesis. To validate these results, we tested nine small molecule inhibitors of five of those putative targets and found a dose-dependent effect on four of them, notably those targeting HSP90 and PI3K. These results highlight the predictive power of our personalised Boolean models and illustrate how they can be used for precision oncology.

WebMaBoSS: A Web Interface for Simulating Boolean Models Stochastically.

WebMaBoSS is an easy-to-use web interface for conversion, storage, simulation and analysis of Boolean models that allows to get insight from these models without any specific knowledge of modeling or coding. It relies on an existing software, MaBoSS, which simulates Boolean models using a stochastic approach: it applies continuous time Markov processes over the Boolean network. It was initially built to fill the gap between Boolean and continuous formalisms, i.e., providing semi-quantitative results using a simple representation with a minimum number of parameters to fit. The goal of WebMaBoSS is to simplify the use and the analysis of Boolean models coping with two main issues: 1) the simulation of Boolean models of intracellular processes with MaBoSS, or any modeling tool, may appear as non-intuitive for non-experts; 2) the simulation of already-published models available in current model databases (e.g., Cell Collective, BioModels) may require some extra steps to ensure compatibility with modeling tools such as MaBoSS. With WebMaBoSS, new models can be created or imported directly from existing databases. They can then be simulated, modified and stored in personal folders. Model simulations are performed easily, results visualized interactively, and figures can be exported in a preferred format. Extensive model analyses such as mutant screening or parameter sensitivity can also be performed. For all these tasks, results are stored and can be subsequently filtered to look for specific outputs. This web interface can be accessed at the address: https://maboss.curie.fr/webmaboss/ and deployed locally using docker. This application is open-source under LGPL license, and available at https://github.com/sysbio-curie/WebMaBoSS.

A high-risk retinoblastoma subtype with stemness features, dedifferentiated cone states and neuronal/ganglion cell gene expression.

Retinoblastoma is the most frequent intraocular malignancy in children, originating from a maturing cone precursor in the developing retina. Little is known on the molecular basis underlying the biological and clinical behavior of this cancer. Here, using multi-omics data, we demonstrate the existence of two retinoblastoma subtypes. Subtype 1, of earlier onset, includes most of the heritable forms. It harbors few genetic alterations other than the initiating RB1 inactivation and corresponds to differentiated tumors expressing mature cone markers. By contrast, subtype 2 tumors harbor frequent recurrent genetic alterations including MYCN-amplification. They express markers of less differentiated cone together with neuronal/ganglion cell markers with marked inter- and intra-tumor heterogeneity. The cone dedifferentiation in subtype 2 is associated with stemness features including low immune and interferon response, E2F and MYC/MYCN activation and a higher propensity for metastasis. The recognition of these two subtypes, one maintaining a cone-differentiated state, and the other, more aggressive, associated with cone dedifferentiation and expression of neuronal markers, opens up important biological and clinical perspectives for retinoblastomas.

Network-Based Integration of Multi-Omics Data Identifies the Determinants of miR-491-5p Effects.

The identification of miRNAs' targets and associated regulatory networks might allow the definition of new strategies using drugs whose association mimics a given miRNA's effects. Based on this assumption we devised a multi-omics approach to precisely characterize miRNAs' effects. We combined miR-491-5p target affinity purification, RNA microarray, and mass spectrometry to perform an integrated analysis in ovarian cancer cell lines. We thus constructed an interaction network that highlighted highly connected hubs being either direct or indirect targets of miR-491-5p effects: the already known EGFR and BCL2L1 but also EP300, CTNNB1 and several small-GTPases. By using different combinations of specific inhibitors of these hubs, we could greatly enhance their respective cytotoxicity and mimic the miR-491-5p-induced phenotype. Our methodology thus constitutes an interesting strategy to comprehensively study the effects of a given miRNA. Moreover, we identified targets for which pharmacological inhibitors are already available for a clinical use or in clinical trials. This study might thus enable innovative therapeutic options for ovarian cancer, which remains the leading cause of death from gynecological malignancies in developed countries.

Personalized logical models to investigate cancer response to BRAF treatments in melanomas and colorectal cancers.

The study of response to cancer treatments has benefited greatly from the contribution of different omics data but their interpretation is sometimes difficult. Some mathematical models based on prior biological knowledge of signaling pathways facilitate this interpretation but often require fitting of their parameters using perturbation data. We propose a more qualitative mechanistic approach, based on logical formalism and on the sole mapping and interpretation of omics data, and able to recover differences in sensitivity to gene inhibition without model training. This approach is showcased by the study of BRAF inhibition in patients with melanomas and colorectal cancers who experience significant differences in sensitivity despite similar omics profiles. We first gather information from literature and build a logical model summarizing the regulatory network of the mitogen-activated protein kinase (MAPK) pathway surrounding BRAF, with factors involved in the BRAF inhibition resistance mechanisms. The relevance of this model is verified by automatically assessing that it qualitatively reproduces response or resistance behaviors identified in the literature. Data from over 100 melanoma and colorectal cancer cell lines are then used to validate the model's ability to explain differences in sensitivity. This generic model is transformed into personalized cell line-specific logical models by integrating the omics information of the cell lines as constraints of the model. The use of mutations alone allows personalized models to correlate significantly with experimental sensitivities to BRAF inhibition, both from drug and CRISPR targeting, and even better with the joint use of mutations and RNA, supporting multi-omics mechanistic models. A comparison of these untrained models with learning approaches highlights similarities in interpretation and complementarity depending on the size of the datasets. This parsimonious pipeline, which can easily be extended to other biological questions, makes it possible to explore the mechanistic causes of the response to treatment, on an individualized basis.

Gene- and pathway-level analyses of iCOGS variants highlight novel signaling pathways underlying familial breast cancer susceptibility.

Single-nucleotide polymorphisms (SNPs) in over 180 loci have been associated with breast cancer (BC) through genome-wide association studies involving mostly unselected population-based case-control series. Some of them modify BC risk of women carrying a BRCA1 or BRCA2 (BRCA1/2) mutation and may also explain BC risk variability in BC-prone families with no BRCA1/2 mutation. Here, we assessed the contribution of SNPs of the iCOGS array in GENESIS consisting of BC cases with no BRCA1/2 mutation and a sister with BC, and population controls. Genotyping data were available for 1281 index cases, 731 sisters with BC, 457 unaffected sisters and 1272 controls. In addition to the standard SNP-level analysis using index cases and controls, we performed pedigree-based association tests to capture transmission information in the sibships. We also performed gene- and pathway-level analyses to maximize the power to detect associations with lower-frequency SNPs or those with modest effect sizes. While SNP-level analyses identified 18 loci, gene-level analyses identified 112 genes. Furthermore, 31 Kyoto Encyclopedia of Genes and Genomes and 7 Atlas of Cancer Signaling Network pathways were highlighted (false discovery rate of 5%). Using results from the "index case-control" analysis, we built pathway-derived polygenic risk scores (PRS) and assessed their performance in the population-based CECILE study and in a data set composed of GENESIS-affected sisters and CECILE controls. Although these PRS had poor predictive value in the general population, they performed better than a PRS built using our SNP-level findings, and we found that the joint effect of family history and PRS needs to be considered in risk prediction models.

Modeling Progression of Single Cell Populations Through the Cell Cycle as a Sequence of Switches.

Cell cycle is a biological process underlying the existence and propagation of life in time and space. It has been an object for mathematical modeling for long, with several alternative mechanistic modeling principles suggested, describing in more or less details the known molecular mechanisms. Recently, cell cycle has been investigated at single cell level in snapshots of unsynchronized cell populations, exploiting the new methods for transcriptomic and proteomic molecular profiling. This raises a need for simplified semi-phenomenological cell cycle models, in order to formalize the processes underlying the cell cycle, at a higher abstracted level. Here we suggest a modeling framework, recapitulating the most important properties of the cell cycle as a limit trajectory of a dynamical process characterized by several internal states with switches between them. In the simplest form, this leads to a limit cycle trajectory, composed by linear segments in logarithmic coordinates describing some extensive (depending on system size) cell properties. We prove a theorem connecting the effective embedding dimensionality of the cell cycle trajectory with the number of its linear segments. We also develop a simplified kinetic model with piecewise-constant kinetic rates describing the dynamics of lumps of genes involved in S-phase and G2/M phases. We show how the developed cell cycle models can be applied to analyze the available single cell datasets and simulate certain properties of the observed cell cycle trajectories. Based on our model, we can predict with good accuracy the cell line doubling time from the length of cell cycle trajectory.

Dynamical Boolean Modeling of Immunogenic Cell Death.

As opposed to the standard tolerogenic apoptosis, immunogenic cell death (ICD) constitutes a type of cellular demise that elicits an adaptive immune response. ICD has been characterized in malignant cells following cytotoxic interventions, such as chemotherapy or radiotherapy. Briefly, ICD of cancer cells releases some stress/danger signals that attract and activate dendritic cells (DCs). The latter can then engulf and cross-present tumor antigens to T lymphocytes, thus priming a cancer-specific immunity. This series of reactions works as a positive feedback loop where the antitumor immunity further improves the therapeutic efficacy by targeting cancer cells spared by the cytotoxic agent. However, not all chemotherapeutic drugs currently approved for cancer treatment are able to stimulate bona fide ICD: some commonly used agents, such as cisplatin or 5-fluorouracil, are unable to activate all features of ICD. Therefore, a better characterization of the process could help identify some gene or protein candidates to target pharmacologically and suggest combinations of drugs that would favor/increase antitumor immune response. To this end, we have built a mathematical model of the major cell types that intervene in ICD, namely cancer cells, DCs, CD8+ and CD4+ T cells. Our model not only integrates intracellular mechanisms within each individual cell entity, but also incorporates intercellular communications between them. The resulting cell population model recapitulates key features of the dynamics of ICD after an initial treatment, in particular the time-dependent size of the different cell types. The model is based on a discrete Boolean formalism and is simulated by means of a software tool, UPMaBoSS, which performs stochastic simulations with continuous time, considering the dynamics of the system at the cell population level with appropriate timing of events, and accounting for death and division of each cell type. With this model, the time scales of some of the processes involved in ICD, which are challenging to measure experimentally, have been predicted. In addition, our model analysis led to the identification of actionable targets for boosting ICD-induced antitumor response. All computational analyses and results are compiled in interactive notebooks which cover the presentation of the network structure, model simulations, and parameter sensitivity analyses.

Robust and Scalable Learning of Complex Intrinsic Dataset Geometry via ElPiGraph.

Multidimensional datapoint clouds representing large datasets are frequently characterized by non-trivial low-dimensional geometry and topology which can be recovered by unsupervised machine learning approaches, in particular, by principal graphs. Principal graphs approximate the multivariate data by a graph injected into the data space with some constraints imposed on the node mapping. Here we present ElPiGraph, a scalable and robust method for constructing principal graphs. ElPiGraph exploits and further develops the concept of elastic energy, the topological graph grammar approach, and a gradient descent-like optimization of the graph topology. The method is able to withstand high levels of noise and is capable of approximating data point clouds via principal graph ensembles. This strategy can be used to estimate the statistical significance of complex data features and to summarize them into a single consensus principal graph. ElPiGraph deals efficiently with large datasets in various fields such as biology, where it can be used for example with single-cell transcriptomic or epigenomic datasets to infer gene expression dynamics and recover differentiation landscapes.

Modelling of Immune Checkpoint Network Explains Synergistic Effects of Combined Immune Checkpoint Inhibitor Therapy and the Impact of Cytokines in Patient Response.

After the success of the new generation of immune therapies, immune checkpoint receptors have become one important center of attention of molecular oncologists. The initial success and hopes of anti-programmed cell death protein 1 (anti-PD1) and anti-cytotoxic T-lymphocyte-associated protein 4 (anti-CTLA4) therapies have shown some limitations since a majority of patients have continued to show resistance. Other immune checkpoints have raised some interest and are under investigation, such as T cell immunoglobulin and ITIM (immunoreceptor tyrosine-based inhibition motif) domain (TIGIT), inducible T-cell costimulator (ICOS), and T cell immunoglobulin and mucin domain-containing protein 3 (TIM3), which appear as promising targets for immunotherapy. To explore their role and study possible synergetic effects of these different checkpoints, we have built a model of T cell receptor (TCR) regulation including not only PD1 and CTLA4, but also other well studied checkpoints (TIGIT, TIM3, lymphocyte activation gene 3 (LAG3), cluster of differentiation 226 (CD226), ICOS, and tumour necrosis factor receptors (TNFRs)) and simulated different aspects of T cell biology. Our model shows good correspondence with observations from available experimental studies of anti-PD1 and anti-CTLA4 therapies and suggest efficient combinations of immune checkpoint inhibitors (ICI). Among the possible candidates, TIGIT appears to be the most promising drug target in our model. The model predicts that signal transducer and activator of transcription 1 (STAT1)/STAT4-dependent pathways, activated by cytokines such as interleukin 12 (IL12) and interferon gamma (IFNG), could improve the effect of ICI therapy via upregulation of Tbet, suggesting that the effect of the cytokines related to STAT3/STAT1 activity is dependent on the balance between STAT1 and STAT3 downstream signalling.

Trajectories, bifurcations, and pseudo-time in large clinical datasets: applications to myocardial infarction and diabetes data.

BACKGROUND: Large observational clinical datasets are becoming increasingly available for mining associations between various disease traits and administered therapy. These datasets can be considered as representations of the landscape of all possible disease conditions, in which a concrete disease state develops through stereotypical routes, characterized by "points of no return" and "final states" (such as lethal or recovery states). Extracting this information directly from the data remains challenging, especially in the case of synchronic (with a short-term follow-up) observations. RESULTS: Here we suggest a semi-supervised methodology for the analysis of large clinical datasets, characterized by mixed data types and missing values, through modeling the geometrical data structure as a bouquet of bifurcating clinical trajectories. The methodology is based on application of elastic principal graphs, which can address simultaneously the tasks of dimensionality reduction, data visualization, clustering, feature selection, and quantifying the geodesic distances (pseudo-time) in partially ordered sequences of observations. The methodology allows a patient to be positioned on a particular clinical trajectory (pathological scenario) and the degree of progression along it to be characterized with a qualitative estimate of the uncertainty of the prognosis. We developed a tool ClinTrajan for clinical trajectory analysis implemented in the Python programming language. We test the methodology in 2 large publicly available datasets: myocardial infarction complications and readmission of diabetic patients data. CONCLUSIONS: Our pseudo-time quantification-based approach makes it possible to apply the methods developed for dynamical disease phenotyping and illness trajectory analysis (diachronic data analysis) to synchronic observational data.